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SUMMARY 


A modified stepwise regression is applied to flight data from a light research 
airplane operating at high angles of attack. The well-known phenomenon referred to 
as "bucking" or "porpoising" is analyzed and modeled using both power series and 
spline expansions of the aerodynamic force and moment coefficients associated with 
the longitudinal equations of motion. The resulting models are validated by numer- 
ically integrating them using initial flight conditions and flight control inputs. 
In addition, a one -degree -of -freedom Van der Pol model is used to help explain the 
oscillatory behavior, and the possible existence of hysteresis in the lift curve is 
demonstrated. 


INTRODUCTION 

An analysis of airplane flight data which exhibits seemingly spontaneous short- 
period longitudinal oscillations occurring at high values of is presented. Such 

behavior has popularly been referred to as "bucking" or "porpoising." Phillips 
(ref. 1) developed several possible models to simulate such behavior on the hybrid 
computer. His models involved hysteresis loops in the aerodynamic force and moment 
coefficients in the region of maximum lift. Other authors (refs. 2 and 3) have 
developed models for self -excited longitudinal motion in the deep-stall condition. 
This paper addresses another example of bucking or porpoising behavior. Oscillations 
were encountered at high values of C L by a light airplane modified specifically for 
high a operation. A published report (ref. 4) on wind-tunnel tests of a model of a 
similar (but smaller) airplane with the same modification indicates that maximum lift 
is at a significantly higher angle of attack. By applying a stepwise regression 
technique to flight data from a nonlinear operating regime of a light single-engine 
research airplane, mathematical models are synthesized and aerodynamic parameters are 
evaluated. Two methods, one utilizing partitioning or binning of the data, and the 
other employing spline basis functions, are used in the model structure determination 
and parameter estimation process. A brief description of each method is contained 
herein. In addition, the Van der Pol equation is used to help illustrate and explain 
the oscillatory motions of the airplane. The flight data are also analyzed in an 
effort to detect possible hysteresis in the lift curve. 

A better understanding of nonlinear aerodynamic phenomena and the analysis of 
such phenomena are important in aviation safety and in the synthesis of flight con- 
trol laws for nonlinear operating regimes (ref. 5) . 


SYMBOLS 


a X ' a y' a Z acceleration along longitudinal, lateral, and vertical body axes, 
respectively, g units 


b 


span, m 


C L lift coefficient, L/qS 

C. maximum lift value of lift coefficient, L/qS 
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rolling-moment coefficient, M x /qSb 
pitching -moment coefficient, My/qSc 
defined in appendix 
yawing-moment coefficient, M z /qSb 
longitudinal-force coefficient, F x /qS 
side-force coefficient, Fy/qS 
vertical-force coefficient, F z /qS 
center of gravity 
mean aerodynamic chord, m 
F-statistic 

F-statistic used in partial F-test 

force along longitudinal, lateral, and vertical body axes, respectively, N 
acceleration due to gravity, m/sec^ 

moment of inertia about longitudinal, lateral, and vertical body axes, 
respectively, kg-m 2 

product of inertia, kg-m^ 

cost function 

number of spline knots 

lift, N 

rolling, pitching, and yawing moments, respectively, N-m 
mass , Kg 

number of data points 

body axis roll rate, rad/sec or deg/sec 

body axis pitch rate, rad/sec or deg/sec 
l 2 

= ~ pV , dynamic pressure. Pa 

squared nultiple correlation coefficient 

body axis yaw rate, rad/sec or deg/sec 
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wing area, ra 2 

variance of parameter estimate, 0^ 
time, sec 

co«qc>onent of velocity along longitudinal body axis, m/sec 
total airspeed, m/sec 

component of velocity along lateral body axis, m/sec 

component of velocity along vertical body axis, m/s ec 

jth element of n x 1 vector of independent variables 

n x 1 vector of dependent variables 

dummy variable 

angle of attack, rad or deg 

value of angle of attack corresponding to ith spline knot, rad or deg 

(a - a^) m (a > o^) 

0 (a < a^) 

angle of sideslip, rad or deg 

= a - a Q , rad or deg 

difference between trim angle of attack and angle of attack at Which pitch 
damping changes sign 

aileron deflection, rad or deg 

stabilator deflection (positive for trailing edge down), rad or sec 
rudder deflection, rad or deg 

damping coefficient in Van der Pol differential equation 

pitch angle, rad or deg 

(n + 1) x i vector of unknown parameters 

jth element of vector of unknown parameters 

air density, kg/m^ 

bank angle, rad or deg 

characteristic oscillation frequency, rad/sec 
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Subscript: 


trim value 


Superscripts: 


derivative with respect to time 


optimal estimate 


Derivatives : 


q 5q c/2V 


qa 5q 5a c/2V 


~m 5a 
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• 5a c/2V 
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FLIGHT VEHICLE AMD DATA ACQUISITION 

The test-flight vehicle was a single-engine, low-wing, four-seat, light air- 
plane. Flight control was provided by the throttle, stabilator (all movable tail), 
aileron3, and rudder. Flaps were available, but were not used in the tests consid- 
ered in this report. The wing leading edge had been modified by extending the out- 
board half of each leading edge as indicated in figure 1 and in modification R for 
the model of a two-place airplane in reference 4. This was one of several leading- 
edge modifications that were made on this airplane to research the stall/spin charac- 
teristics of such a modified airplane. The flight data analyzed herein were obtained 
during standard stability and control parameter estimation flights. In these flights, 
the pilot attempts to trim the airplane to some steady-state reference flight condi 
tion from which he can initiate some perturbation by the movement of controls# All 
data analyzed in this report are assumed to have equilibrium initial conditions 
(p=q*r*a=|5=p = v = w=*u = 0). Data were measured by rate gyros for roll, 
pitch, and yaw rates; roll and pitch angles were measured by attitude gyros. Angles 
of attack and si * slip were measured by wind vanes mounted on booms on each wing tip 
as shown in figure 1 . Linear accelerations were measured by accelerometers located 
close t the airplane center of gravity with three mutually orthogonal axes pointing 
in the directions of the longitudinal, lateral, and vertical body reference axes of 
the airplane. Control displacements were measured by potentiometers located close to 
the respective control surface so as to eliminate time delays and inaccuracies in the 
measurements of displacement due to control-cable stretch. 

Data were recorded onboard the test airplane as analog voltage signals. The 
records of these voltages were filtered by a 6-Hz low-pass filter, and digitized to a 
sample rate of 20 per second after each flight. The digitized 20-sample-per-second 
data were then corrected for c.g. offset of the instruments, upwash and bias error on 
the wind vanes, and bias errors in the accelerometer and the angular-rate gyros. 

ANALYTICAL PROCEDURE 
System Identification 

The system identification problem was defined by Zadeh (ref. 6) as "determina- 
tion, on the basis of observation of input and output, of a system within a speci- 
fied class of systems to which the system under test is equivalent." For the system 
identification of an airplane operating at low angle of attack, the mathematical 
model structure of aerodynamic forces and moments is linear. Hence, the identifica- 
tion problem reduces to a parameter estimation problem. However, at high angle of 
attack and in near-stall operating regimes, the form of the aerodynamic forces and 
moments must be determined before estimating corresponding parameter values. The 
general form of the mathematical model structure for the aerodynamic force and moment 
coefficients can be written as 


y(t) » 0 n + 8.x (t) + 0 x (t) + ••• 0 x (t) (1) 

oil 22 n n 

where 

y(t) aerodynamic force or moment coefficient (C ,C ,C ,C ,C ,C ) at time t 

X Y Z m x n 

x.i(t) airplane state plus control variables (a,q,P,p,r,6 e ,6 a ,6 r ) and their com- 
binations at time t (j - 1, 2, ..., n) 
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0^ airplane stability and control coefficients (j * 1, 2, ..., n) 

Qq constant reflecting any initial steady-state condition 

Thus, although it is assumed for the identification problem at low angles of attack 
that only first-order terms in the state and control variables are required for the 
m^del^the identification problem at high angles of attack can require terms such as 
a , a , qa, and 6 a* Hence, part of the identification problem is the proper 
choice of these higher-order terms for which parameters must be estimated. For exam- 
ple, consider the vertical-force coefficient C 2 * For the linear model, 


c z = 


+ C Ao + C 


£C 

Z 2V 


+ C_ 


A6 


( 2 ) 


where Act = a - a and A6^ = 6_ - 6„ for steady state (v = p„ = at = r = * = 0) 

o e e e Q o o td o o 

initial flight conditions. A corresponding nonlinear model for high-angle -of -attack 
maneuvers might be written as 


c z - 


+ C Ao + C 


S£ + c 

Z 2V Z. 

q 6 


46 e * C z 2 

a 


(Aa) + C„ 


Aa 


qc 

2V 


(3) 


The specified class of models from which a nonlinear model may be chosen is given by 
the candidate model variables (table 1). Table 1 is a set of influence variables for 
longitudinal motion. The first column is simply the linear model variables a, q, 
and the linear control variable 6 0 , on which the short-period longitudinal motion 
depends. The second column is the variation of the linear terms with changing angle 
of attack. The terms an( j p2 a i n the third column allow for aerodynamic coupl- 

ing to lateral motion, and the last column contains the terms allowing strong non- 
linear variation with angle of attack, the main longitudinal independent variable. 

The list of variables in table 1 can be expanded or changed to investigate any parti- 
cular behavior that is a function of a, q, and 6 e » 

Once a class of candidate model variables is specified, the problem still 
remains to select the proper subset of that class. For the work described in this 
report, a semiautomated procedure, the stepwise regression, is employed. The step- 
wise regression parameter estimates are actually linear regression parameter esti- 
mates predicated on the minimization of the equation error, whose mean square is 
given by the cost function 


J = i E - y(i)i 2 u> 

i = 1 


where 

N 

y(i) * e 0 + E x (i) e 

j = 1 J 
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and the values of 0. are the equation error optimal parameter estimates (j ■ 0, 
1# . n) . In the Stepwise regression, the x. 's are added to the regression 
according to their effectiveness in reducing tne equation error cost function J. 
This value is represented by the partial F value F of the coefficient as 
follows: P 


F 

P 


0 


I 


s 2 (0 j ) 


( 5 ) 


A A a A 

where 0. is the estimate of 0* and s^(0. ) is the variance of 04 . The value of 
the stepwise regression in system identification is that it not only selects a model 
from a large class of candidate model variables but it simultaneously estimates the 
corresponding parameters for that model. 

When applying a stepwise regression to flight data, it should be remembered that 
the data are corrupted by measurement noise, and perhaps by process noise, which 
appears in the form of modeling error. Hence, the regression, with fewer independent 
variables than data measurements, always adds additional terms in order to reduce J. 
Though all these terms tend to fit the data better, some terms simply fit the noise. 
Hence, some criterion is needed to determine the point at which the adding of vari- 
ables should cease. Experience has shown that several criteria should be employed to 
ensure an adequate model that fits the data, has good prediction capabilities, and 
makes sense physically. (See ref. 7.) These criteria include; 

1. Fp should be greater than 5. 

2. F should be a maximum. 

3. R should be close to 100 percent. 

4. The residual sequence should be statistically equivalent to white noise. 

This determination is made by observing the autocorrelation function for 
the residual sequence. 

No satisfactory single performance index has been written as a function of these 
separate criteria? hence, the selection of an adequate model is still somewhat 
subjective. 


Partitioning of Data 

It was demonstrated in reference 8 that for large-amplitude longitudinal maneu- 
vers the partitioning of data as a function of angle of attack provided a finer 
resolution of stability and control parameters for the maneuver range than could be 
achieved by analyzing the entire data string intact. This partitioning, or binning, 
technique can be applied as follows. Consider a data set consisting of the measure- 
ments of several variables of 400 time points each. The time history of one such 
variable, a, is shown in figure 2. Several bins can be created, each containing 
data corresponding to a given range of a. The lines drawn parallel to the abscissa 
in figure 2 represent the bin boundaries. For example, all data corresponding to 
0° < a < 4° would be associated with the first bin, cu v .a for 4° < a < 8° would be 
associated with the second bin, and so forth, as shown in figure 2. Then each bin. 
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which now contains data corresponding to a limited range of angle of attack, can be 
analyzed for aerodynamic stability and control derivatives corresponding to that 
angle-of -attack region. 

Hiere are, however, several problems associated with identifying parameters from 
binned da La. The first is that by arbitrarily assigning bin widths, the a region 
may not be d. vided finely enough. However, one might also divide the region so that 
there are too little data in each bin for a proper analysis. Moreover, the location 
of the hin enopoints might obscure an important feature of the data if the endpoints 
were shifted a degree. Hence, a second technique was applied to the data in an 
attempt to alleviate these potential problems. This second technique allows the data 
to "bin itself" through the use of spline functions. This approach is more general 
i*i applicability than the simple binning and can be described as follows. The range 
of the independent variable which is most important in the determination of the 
dependent variable is partitioned into several subsets, each having support on less 
of the range than the previous subset. For example, the force coefficient is 

mainly dependent on a. Hence, if a * {z|a < z < b>, then the a range, [a,b], is 
divided according to the spline basis functions as follows: 


(a - 



f (a - a L ) m (a > c^) 

|^0 (a < a^) 


The values of are called knots. An example of the "+" function is given in 

figure 3. The four knots in this figure are at a = 2°, 4°, 6°, and 8°. Hence, 

(a - ® = 1 for a > =2°, and (a - a 1 = 0 for a < . Similarly, 

(a - a 2 ) + « 1 for a > 4°, and (a - = 0 for a 4°, and so forth, for the 

rest of the "+" functions. If the order of the "+" function, denoted by the super- 

script m is other than zero, say 2, then (a - * (a - otj) for a > a-j 

and ( a - OLj ) l - 0 for a < • For the analysis in this paper, the knots were 


placed every 0.5° between a *= 10° and a = 18°. By placing knots every 0.5° and 
then allowing the stepwise regression to determine significant terms, the problem 
of choosing too narrow a bin in the binning technique J , avoided. The longitudinal 
force and moment coefficients were then written as follows: 


c x = 


K 

C X,o + °X « + 

a i=l 


'X (a ‘ a iK + C X M 
“i q 


K 

£ 

i=1 


: (a-a.)°3£- 

X i + 2V 

q i 



+ 



V a " 

7 


«7>J 



6 e (a - 


13 


„ v0 

a 13 + 
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'2 “ C z,o + C z a + £ C z (« - o i ) + ♦ C z % + r C z (a 

a i»1 a. q i-t q. 

1 t 

* C Z ^ 6 . ♦ - « 7 »i ( C * 6< ] 13 6 * (0 ‘ “' 3 1 ' 


»! 


> ■ C ».» * V ‘ ‘ “3’* * % » * £ \‘* 

* V* * (%.)/•“ - *** * (S.) 13 Va ’ “’ 3,? 



3c 

2V 



as. 

2V 


Only two k a, a and a at a - 13* and a - 16°* were chosen for the ele- 
vator because of a 7 lack of computer memory, the locations of these two knots were 
chosen to allow for what appeared to be the three ma^or areas of elevator 
effectiveness. 

After a model has been selected and parameters have been estimated, the predic- 
tive capabilities of the model must be tested, this is known as validation of the 
model, the method of validation used in this report is the numerical integration of 
the equations of motion, this method uses the estimated model and a fourth-order 
Runge-Kutta integration scheme. Starting with the initial conditions of the flight 
for which a given model was estimated, the longitudinal output variables a, q, 0, 
and V are confuted and compared with the actual flight values for these variables. 
The robustness of the model can be validated by applying it to another flight with 
only similar initial conditions. 


RESULTS AND DISCUSSION 
Flight Data 

A typical "flight," or ejqperiment, for the determination of stability and con- 
trol parameters consists of several "runs." Each run consists of about 3 minutes of 
data commencing when the airplane achieves some predetermined altitude (usually about 
2000 m for a light single-engine airplane). For idle power experiments, after 
attaining the predetermined altitude, the throttle is closed. The airplane is then 
trimmed to a desired angle of attack. A control input perturbs the airplane from the 
trim condition with the perturbation damping out in about 4 seconds for short-period 
longitudinal disturbances. After the disturbance has damped, the airplane is trimmed 
again to some desired angle of attack and is perturbed again by a control input. 

In this manner, five to seven perturbations can be carried out during a single run. 
Each run within a flight is designated by an integer - successive runs are designated 
by successive integers. Each perturbation within a run is designated by a decimal 
point and digit appended on the run number. For example, run 7.1 indicates the first 
perturbation in the seventh run of a given flight; run 7.2 is the second perturbation 
executed in the seventh run. A run ends when the airplane approaches some predeter- 
mined minimum altitude (about 1000 m for the experiment reported herein). The air- 
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plane then climbs back to the initial altitude (taking about ID minutes) and proceeds 
with the next run* The attitude gyros are caged (locked) between runs and uncaged at 
the commencement of a run* 


Ficr: * presents traces for the longitudinal forces, certain lateral-motion 
variahH , an "he longitudinal control variable. The figure shows that after being 
pertu' -d from trim angle of attack of about 13.6°, the stabilator is brought back 
to a position th*iC is displaced by 1° from its trim position* While the stabilator 
remains this position, the induced longitudinal oscillation does not damp out, but 
continues 3‘; constant amplitude until the stabilator is returned to the -5° initial 
condition, .he time histories of the lateral-motion variables indicate that coupling 
between the lateral and longitudinal modes is not involved* 

Figure 5 presents traces of the longitudinal-motion variables and control dis- 
placements for runs 7,1 through 7*6 and runs 8.1 through 8.5, all of which exhibit 
similar undamped oscillatory behavior. 

Confirmation of the oscillatory region is shown in the plot of run 11 in fig- 
ure 6. In this run, the airplane was trimmed to an angle of attack below that at 
which the undamped oscillations occurred. The stabilator deflection was increased 
in increments of less than 1°. Hie oscillation is centered around a = 14° and 



Stabilator deflection to trim is plotted against trim angle of attack in fig- 
ure 7. Hiis figure can be used to check the actual stability of the trim conditions 
in the 14° to 16° angle-of -attack range. From figure 7, a loss of stabilator effec- 
tiveness to trim is indicated between a * 14° and a * 17.5°. 


For each run in figure 7, an uncertainty range for 6 e was calculated by aver- 
aging the 10 sample points for stabilator deflection before and the 10 sample points 
for stabilator deflection after the pilot noted that the airplane was trimmed. In 
the region above a =» 14° and below a - 18°, there is a very large uncertainty in 
the value of stabilator to trim. This uncertaincy is shown by the bars about the 
stabilator to trim values in this region. 

Hie lower bound on this region is further defined by runs in which the airplane 
trimmed at angles of attack near 14° but pitched nose down instead of nose up. 

Hie oscillations induced in this manner do damp out. Hence, one is left with a 
steady-state short-period longitudinal oscillation which is easily induced from trim 
angles of attack in a region bounded by a = 14° and a = 17° and with stabilator 
deflections between -7° and -9° with the throttle set at idle power. 


Comparison of Binned and Spline Analyses 

Figure 8 presents the results of the stepwise regression analysis applied to 
three of the oscillatory flights. Hie circles represent linear parameters estimated 
from the binned analysis. Bach circle is plotted according to the mean value of a 
for data in the bin it represents. Hie short-line segments in figure 8 represent the 
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a lope of the parameters. These slopes can be identified when ‘ second-order term 
such as 


C 


Z 




a 



1 J&_ 

2 da 



( 7 ) 


enters the regression equation. 

As with the values of the parameters themselves# each slope for a para^neter is 
plotted at the mean value of a for the bin in which it is identified. The ordinate 
value at which a slope is plotted is determined by extrapolating a line with the 
identified slope from the identified parameter value at the airplane trim angle of 
attack. (See fig. 8(a)# for mple.) On the plot of C_ as a function of a 

Z q 

in this figure# a circle is plotted at a * 11.3° and C_ = -22. However# for the 

q 

next bin# with a mean angle of attack of 12.6°# a slope is plotted* This slope 

(-390 rad" 1 ) is given by the value of C z # which is the coefficient of the a ^ 

ga 2V 

term identified for this bin. The value of C„ for the data in this bin is -47 but 

Z q 

is referred to a Q for the run which is 16.4°. Hence# the extrapolation equation 
for C 2 to the bin which has a mean value of a = 12.6° is 


r_ vu> = -47 - 390 (a - a n ) (8) 

Z q ° 


where a * 0.2862 radians and a is in radians. Experience has shown that such a 
nonzero 9alue for a slope of a derivative indicates a uon of change in the deriva- 
tive value with respect to angle of a tack. 

Also plotted in figure 8 are the results from the application of the stepwise 
regression with the spline candidate model. The dat/ are not £ irtitioned prior to 
the application of the spline model. Any change in parameter values with angle of 
attack is noted by the entry of a knot near the region of change. Here again# the 
plot of C„ as a function of a in figure 8(a) provides an example. The "stair- 

Z q 

case" or step-shaped line is the result of a spline analysis in which the stepwise 
regression selected zero-order splines in (a - a^)^ ^ with knots at a = 13.5° 
and 14.5° or a model for C z given (for the entire maneuver) by 

q 

C„ (a) - -23. - 5.5(o - 0.2356)° - 5.0(r: - 0.2530)° (9) 

Z + + 

q 
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where a is in radians. An examination of the plots for each derivative in figure 8 
is included xn the remainder of this section* 

Both the binning and spline technique show C* (figs* 8(a), 8(c), and 8(e)) 

to consist of at least two levels* For all three runs, C„ is about -2*2 between 

Z a 

a * 12° and a m 15.5°. Above a * 15*5°, C_ decreases in absolute value to 

Z a 

about -1. this corresponds to a flattening of the lift curve at about a * 15*5 P . 

In addition run 7*6, which contains more data from the low-angle-of -attack region, 

gives rise to a third level* At this level, C » -4*5 for a < 12°* Overall, 

z 

a 

C 7 is strongly and consistently identified* 
a 

the derivative C 2 (figs* 8(a), 8(c), and 8(e), decreases in the region 

q 

between a * 12° and a * 16° for all three runs* Run 7*6 contains an extra step 
in the low-angle -of -attack range* This is probably caused by a relatively large 
amount of data in that region, as was the case with the parameter C„ • 

the derivative C* (figs. 8(b), 8(d), and 8(f)) is the third strongly and 

“q 

consistently identified derivative* Both the binned and the spline analyses indicate 
positive values of C* between a * 12° and a * 17°* Run 7.6 (f5g* 8(f)) again 

q 

gives the best results in the low-angle -of -at tack range, where the spline technique 

showed that O' = -13 for a < 10*5°-* 

01 

q 

The derivative C* is not identified consistently across the three runs* In 
m a 

run 7.1 , the binned result is scattered, and the spline picks out no structure* The 
fact tint the spline result simply averages the binned results indicates that the 
scatter is probably not a function of the choice of bin width and/or bin boundaries* 

Run 7*3 yields a bi-level value of C* , and the binned-value analysis indicates a 

m a 

drop in C’ from -1.7 to -2.8 at a * 16*6°. The spline analysis of run 7.3 shows 


a similar drop (but to -3.6) at a * 18°* The difference in the angle of attack at 
which the C f change occurs could be explained by aliasing due to bin width and bin 
m a 

boundaries* Run 7,6 shews another slightly different result. The data for the low- 
angle-of -attack range have again resulted in better identification of at 

a < 12°* The binned analysis indicates a trend from about C* * -1 at a » 10° 

% 

to C' «* -2 for a > 16°. The trend is shown to support a type of double break in 
m 

a 

when analyzed by the spline technique. Bach of the breaks identified by the 

CL 

spline model gives C* - - 1. 

% 

Hie derivative C„ is estimated consistently for a > 16° by both techniques 

% 

as C_ = -1,2* At lower angles of attack and in the re = 1 2° to a = 16° region. 
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the three runs yield different values* Run 7*6 again yields the result aost in 
agreement with other data from the low-angle-of -attack region. An aliasing problem 
could exist with the estimation of C„ since, as pointed out previously, only two 

\ 

knots (at a * 1 3° and a = 16°) could be chosen for 6 e * 


The derivative C* showed the same characteristics as C_ • Run 7.6 pre- 

”6 e X 

sented the best low-angle-of -attack values, and all three runs yielded similar values 
of **-1.2 for C 1 for a > 16°. 

"6 e 


Interpretation of Derivatives 


Because each derivative can be interpreted as a coefficient in the differential 
equations describing the motion of the airplane, each derivative can be interpreted 
as representing some characteristic behavior in that motion. The derivative C* 

ra a 

represents the static margin of the airplane. A large negative increase, as seen in 
the analysis of runs 7.3 and 7.6, usually indicates wing stall and the resultant 
shift of tlie aerodynamic center. A positive value of C' , which was identified in 

m q 

all three runs for the region between a * 12° and a » 17°, indicates positive 
pitch damping. Positive pitch damping indicates an unstable or divergent pitch rate 


in this angle of attack region. The decreasing value of 


can be interpreted 


as indicating a flattening of the lift-curve slope. This can also be interpreted as 
a loss in damping in a* The derivatives C,- and C* indicate control effec- 


tiveness. At low angles of attack, as indicated from run 7.6, the ratio of C! 

me 

C„ should be proportional to the moment arm for the stabilator. However, for 


to 


a > 10°, the values of C' 


and C n 


approach one another. These values are 


% % 

predicted on only a few data points overall, particularly in the oscillatory region, 

where the stabilator was held at 6 for about 90 percent of the time. 

e,o 


Validation 

After a model has been selected by an identification technique, it still must be 
validated in some manner. Several methods of validation are commonly used in system 
identification. One method compares the model derived by using one technique with 
th* derived by using another technique. For example, the model selected by using 
an equation error technique can be compared with one derived by using a maximum- 
l'kellhood output error approach. Another method of validation is the testing of 
the predictive capabilities of the model. Given the initial conditions and forcing- 
function time history, the model equations can be numerically integrated, and the 
output variables can be compared with the actual measurement time histories of the 
system output variables. Another method of validation, and one which demonstrates 
the robustness of the model, is the combination of the initial conditions and 
forcing-function time histories of an experiment similar to (but not the same as) 
that used for the identification. Again, the output of the model should be close to 
the output of the actual system. 
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The second method, that of predictive ability, was used to validate the Models 
derived fro« oscillatory data* this technique is particularly suited to the problea 
considered in this report, because the Model was determined by minimizing the Mean* 
square error between the measured and computed acceleration* The integrated equa- 
tions of motion, however, compare output variables - total velocity V, angle of 
attach a, pitch rate q, and pitch angle 6* An inadequate model would fit the 
data to the acceleration time history, but would not, from initial conditions and a 
forcing-function time history, correctly predict the motion of the system* 

The flight time histories for selected longitudinal variables from run 7*1, as 
well as the simulated time history for that run from the integrated longitudinal 
equations of motion with a linear aerodynamic model, are plotted in figure 9(a)* 
Figure 9(b) shows simulated results using the spline model estimated from run 7*1* 
This model gives slightly better results for the angle-of -attack prediction, although 
it is no better than the linear model in predicting the time histories of the other 
three variables* Both models fail to match the total airspeed by about 10 percent at 
the worst points* However, the short-period oscillation in the airspeed is predicted 
by both the linear and nonlinear model* 

As a measure of robustness of the nonlinear model, the initial conditions and 
forcing-function time history of run 7.3 were integrated using the nonlinear model of 
run 7*1* The results are presented in figure 9(c), where the simulation result is 
shown with the actual time history of run 7*3* Run 7.3 is a good test of the model, 
as it was excited with a smaller and less persistent stabilator Input, and consists 
only of the free -oscillation response of the airplane* The model seems to be overly 
dependent on the forcing function, in that it predicts the first oscillation well but 
then allows the oscillation to damp too quickly. The model does, however, simulate 
the free oscillation qualitatively. 


Correlations With Theoretical Models 

The comparison of models estimated from flight data with those identified from 
wind-tunnel experiments or with those synthesized theoretically is complicated by 
several factors. These factors include a difference in Reynolds number, possible 
gust effects in flight, less repeatability of exact initial conditions in flight than 
in the wind tunnel, and, very importantly, an inability to separate the plunging 
motion of the airplane from its pitching motion in flight* The inability to separate 
these two motions is demonstrated in figure 10, where both 4 (plunging) and q are 
plotted against time for one of the oscillatory flights* The values of 4 for fig- 
ure 10 were calculated by applying a cubic spline under tension to the angle-of- 
attack time history, which was indicated by the wind-vane angle-of -attack sensor on 
the airplane. The correlation between these two variables means that aerodynamic 
derivatives with respect to a and q cannot be estimated separately, but rather 
some combination must be estimated. The actual combinations that are estimated are 
derived for the pitching-moment derivatives in the appendix. Within the framework of 
the above difficulties, the models estimated from the flight data were compared with 
two theoretical models* The first type of model was developed by Phillips (ref. 1). 
The model demonstrates that the dynamic effect of 4 near C L mAX of the wing and 
the associated hysteresis in the lift curve can produce an oscillatory longitudinal 
motion* Hence, this first model deals principally with the effect of the wing* The 
second model considered here places the cause of the oscillation mainly on the tail. 
In this second model, when C > 0 in a small a region, a total loss of damping 

m q 

in pitch is Indicated* Hence, a longitudinal oscillation develops to a limit cycle. 
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Both models are considered individually* The Phillips models are of interest in 
light of the results of Winkelmann (ref* 9) and the new wind-tunnel results of 
Newsom, Satran, and Johnson (ref* 4)* ’These results indicate either hysteresis or 
uncertainty in the lift-curve slope in certain angle-of -attack regions* This section 
of the present report is an investigation of existence of hysteretic behavior in the 
actual flight data* 


Phillips aodels* - To apply the stepwise regression for the purpose of estimating 
a model structure similar to the Phillips models, the set of candidate model vari- 
ables was modified to include the variables jfa and \|a a* The l[a and a 
terms were restricted to a > 0, in that it was only for these positive vales of a 
that the augmentation of the lift curve is postulated* The values of a < 0 are 
ignored by the program when creating the data string for the and a terms* 

Of the 11 oscillatory runs that were analyzed, only 3 yielded a parameter pro- 

portional to >[a. In each case, the identified parameter was C „ and C 7 ^ 


a a 


was zero. A typical value for ~ was -12*5 ± 2.0. This value of C z _ 

ia a fa a 

implies an increase in max as a linear function of the difference between flight 
reference angle of attack and stall angle of attack. 


A second feature of the Phillips models was hysteresis in the lift curve. To 
identify possible hysteresis from the flight data, 2 bins were created for each of 
the 11 runs. All data corresponding to a > 0 were put in the first bin. The other 
bin contained data corresponding to a < 0* The stepwise regression was applied to 
each bin for each of the 11 runs. In 10 of the 11 runs, the bin for a > 0 yielded 
a greater lift-curve slope than the bin for a < 0. The results for lift-curve slope 
are given in figure 11. It is clear from the figure that for all runs except 7*1, 
the lift-curve slope is greater for a > 0 than for a < 0* For run 7.1, the two 
slopes are the same. The dashed lines in figure 11 depict higher-order derivative 

information (e.g., C 2 * C- 3' etc *) extracted by the stepwise regression. When 

"a 

this same technique was applied to data from maneuvers by the same airplane trimmed 
at lower angles of attack (a - 4° to a - 5°), there was no significant difference 
in lift-curve slope. 


Finally, approximate values cf C^, calculated using 


C L - -C z * -C 2 - C z (a - a Q ) + (higher-order term where identified) (10) 

o a 

are plotted in figure 12 for models f m the three representative runs - 7.1, 7.3, 
and 7.6. Again, the differences in lift for & > 0 compared with a < 0 for a in 
the oscillatory region are evident. When the maneuvers initiated from lower angles 
of attack were analyzed using a similar binning as a function of a, no differences 
in lift curve were evident. 

In summary, differences in lift-curve value and lift-curve slope appear in data 
from the longitudinal maneuvers initiated from angles of attack near 16°, but not 
in data from similar maneuvers initiated from lower trim angles of attack (a = 4° 
to 5°). Moreover, in the region in which the differences exist, the lift-curve 
slope and lift-curve value are greater for Increasing angle of attack (a > 0). Such 
a result agrees with the Winkelmann (ref. 9) results for a similar wing in a wind- 
tunnel test. 
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Van der Pol modal .- In reference 2, Thomas and Oollingbourne describe a deep- 
stall limit-cycle condition by a Van der Pol- type differential equation* The 
Van der Pci equation has the general form 


z(t) - ct 1 - z 2 (t)] i(t) t z(t) *« 0 


(ID 


The solution to equation (11) is a liV L cycle because of the change in the sign 
of the damping coefficient e(1 - z^(t)i t occurs whenever z(t) traverses the 
values ±1 . Por application to the airplane problem, the perturbation da, in angle 
of attack from a Q , can be substituted for a(t) in equation (11) to yield 

Aa * 02 Aa + © 3 Aa(Aa) 2 + 0^ Aa (12) 

Equation (12) is in a form amenable To solut on for 0^ , 0^, and 0^ by the step- 

wise regression program with the independent variables Aa and Aou Then the fre- 
quency of oscillation of the system can be estimated by 



Also, (Aa) . the Aa value at which dhe damping of the motion changes sign, is 
given by 


( Aa) cr i t 



Applying the stepwise regression to the data from run 7.3 yields * -7*4, 

©2 * -0.52, and 63 = 183. These three valuer* are typical of the other oscilla- 
tory runs. For these values of 6 -j , § 2 # an<3 ©3* ^ sc * 2 * 7 rad/sec, and 
(Aa) cr ^ t 3 ±3.1°. These values are corroborated by the flight data. The interpreta- 
tion of the Van der Pol results is that the aa term provides no damping and in fact 
tends to increase |a - a Q | to 3.1°, when the coefficient of Aa changes sign and 
provides damping. 


Another derivation of the Van der Pol ooeffic nts comes from extracting 
pitching-moment derivatives. If the pitching-nr .ent equation is written as 


X Y q 


4 pV 2 Sc 


[ c; . 


Aa + 


C* 3£ 
m 2” 

q 


'm 2 

a q 


Aa)' 


££ 

2V 


1 


(13) 
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and if a = q 


for a single degree of freedom, then equation (13) 


Am 


pv 2 s^ ba c pV^Sc 

21 n 2V 21 
Y q Y 


c o 

m 2 

a q 


(Aa)' 



Aa 


( 14 ) 


Substituting values from run 7*1, then w^ sc * 2.7 rad/sec and (Aa) cr . t * ±3.0®, 
which is again in agreement with the flight data. Tb validate the single-degree-of- 
freedom Van der Pol model, the model equation was integrated numerically using a 
fourth-order Runge-Kutta integration scheme. A phase-plane plot of the results of 
such an integration using the Van der Pol model from run 7.3 and realistic initial 
conditions of Aa * 0 and q = -0.16 rad/sec is presented in figure 13(a). In all 
parts of figure 13, the sense of increasing time is clockwise. This model reaches 
its limit cycle trajectory in about two complete cycles. However, in figure 13(b), 
a phase-plane plot of the actual flight data indicates that a uuch more rapid tran- 
sition to the liiL.t cycle trajectory (in about one-fourth of a cycle) should be 
expected* Moreover, (Aa) c ^ it corres P °nds to the maximum deviation in a* In order 
to predict this maximum deviation, the damping should change sign at a smaller value 
of Aa. In summary, the Van der Pol single-degree-of-f reedom equation correctly 
predicts the qualitative behavior of the motion but fails to predict correctly the 
variation in angle of attack. 

For completeness, two additional cases were numerically integrated, and their 
phase-plane trajectories were plotted. In figure 13(c), the Van der Pol model from 
run 7.3 is integrated from initial conditions Aa - 0 and q = -0.45 rad/sec. These 
conditions lead to the same limit cycle trajectory as did the conditions that yielded 
the results shown in figure 13(a). Finally, a linear harmonic oscillator estimated 
from the equation a = a was integrated, and the resulting phase-plane trajectory 
was plotted in figure 13(d). Though it satisfies the range requirement for a, the 
harmonic oscillator model fails to satisfy the range requirement for q and, most 
importantly, fails to have the qualitative property of converging from any given 
initial condition to the same limit cycle trajectory. 


CONCLUSIONS 

Several examples of short period longitudinal oscillatory flight at high values 
of lift coefficient have been analyzed. The data were obtained by a light single- 
engine research ai. rplane with outboard leading-edge modifications for operation at 
high angles of attack. A stepwise regression was applied to the flight data to 
determine aerodynamic force and moment model structure and parameter values. TVo 
separate applications of the stepwise regression were made, each with different forms 
of the candidate models. In one application, the candidate model set was comprised 
of spline basis functions with knots at 0.5° increments in angle of attack. The 
other application consisted of applying the stepwise regression with polynomial 
candidate model terms to the data after that data had been partitioned into bins as a 
function of angle of attack. An analysis of the results, including a comparison with 
theory, indicates the following: 
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1 • lhe change in pitching moment with respect to pitch rate C* is definitely 

q 

positive for the angle-of-attack range over which the oscillations occur, there is a 
flattening of the lift curve in the oscillatory region. Hie values of the change in 
pitching moment with respect to angle of attack C 1 , though not consistent, imply 

“a 

the possibility of a stall in the oscillatory region. 

2. Both techniques (polynomial on a partitioned data set and spline) give 
similar results. 

3. Hie models simulate the oscillatory behavior qualitatively but provide too 
much damping. 

4. Analysis of data partitioned for positive and negative rates of change of 
angle of attack gives strong support for the existence of a hysteresis loop in the 
lift curve for the angle-of-attack range covered by the oscillatory motion. 

5. A one-degree-of-freedom Van der Pol model can simulate the general phase- 
plane behavior of the motion but does not properly limit the range of angle of 
attack. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
November 30, 1 982 
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EQUATIONS OP TOT I ON 

The following equations of motion are given in this appendix* 


u * -qw + rv - g sin 8 + C 

2m x 


• „ PV^S 

v * -ru + pw + g cos 0 sin t + JLm r~ C__ 

2m Y 


v *= -pv + qu + g cos 0 cos ♦ + 


( I - I ' I 2 

Y Z I . XZ, . V . pV Sb _ 

r )*_ c i 


* JZ “ I X \ . I XZ. 2 2. . pV 2 Sc 

q = pr l “V"i + v - p,+ 


* | x ~ M , Ixz f* , P v2sb r 

r = "TX / ^7 (P " qr) ~n7 n 

\ z / z z 


0 ■ q cos ♦ - r sin $ 


♦ ■ p + (q sin $ + r c $) tan 0 


/ 2 2 2 
V * /u + v + w 


^ -1 w 

a * tan — 
u 


• w 
a ■ — 


B - sin’ 1 J 


8 * — 
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C' - c + c § — c 
m x m x m* 4m Z x 

6 a 6 a a 6 a 

e e e 


c , = c + c £S£ c 

Va Va m a 4m Va 


c . , c + c £Sc c 

m.2 m 0 2 m* 4m Z 0 2 

p p a p 


C* . = C . + C c 

mi mi m* 4m Z i 
a a a a 


Assuming steady-state (v Q * p Q * q 0 * r o * * 0) initial flight conditions, the 

longitudinal equations become 


u = -qw - 


g 8 ine+£ aT c x 


• a pV 2 S 

w = qu + g cos 9 + C 

2m 2 


q 


PV 2 Sc , 
2I y m 


9 * q 


and the longitudinal output can be written as 


a„ - — (u + qw + g sin 6) 
X g 


^ • 

a„ * — (w - qu - g cos 0) 
2 g 
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and for the equation error form, the aerodynamic coefficients can be written as 
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TABLE 1.- CANDIDATE MODEL VARIABLES 
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Figure 1.- Three-view drawing of test airplane showing outboard 
leading-edge modifications tnot to scale) on top view. 
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Figure 2.- Data partitioning for a large-amplitude longitudinal maneuver. 
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Figure 3.- Regions of support for spline 
basis functions. 
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(a) Run 7.1 • (b) Run 7.2. 

Figure 5.- Time histories of selected longitudinal variables for osci 

flights . 
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(a) Run 7.1 and linear polynomial model estimated from 

run 7.1 . 

Figure 9.- Measured time histories for selected longitudinal 

variables and those computed from model selected by regression 
program . 
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(b) Run 7*1 and spline model estimated from run 7.1. 
Figure 9*- Continued* 
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Figure 10.- Time histories of pitch rate and rate of change 
of angle of attack for an oscillatory run. 
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Figure 11.- Values of lift-curve slope for oscillatory runs 
indicating possibility of lift-curve hysteresis. 
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(c) Van der Pol model estimated from run 7.3 showing 
attainment of limit cycle. 



«, deg 

(d) Harmonic oscillator estimated from run 7.3. 
Figure 13.- Concluded. 


48 



